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ABSTRACT 

The core-cusp problem remains as one of the unsolved discrepancies between observations and 
theories predicted by the standard paradigm of cold dark matter (CDM) cosmology. To solve this 
problem, we perform A^-body simulations to study the nonlinear response of CDM halos to the variance 
of the gravitational potential induced by gas removal from galaxy centers. In this study, we focus on 
the timescale of the gas ejection, which is strongly correlated with stellar activities, and demonstrate 
that it is one of the key factors in determining the dynamical response of CDM halos. The results 
of simulations show that the power-low index of the mass-density profile of the dark matter halo 
correlated with the timescale of the mass loss, and it is flatter when the mass loss occurs over a short 
time than when it occurs over a long time. However, it is still larger than typical observational values; 
in other words, the central cusp remains for any mass loss model in the simulations. Moreover, for 
the slow mass-loss case, the final density profile of the dark matter halo recovers the universal density 
profiles predicted by the CDM cosmology. Therefore, mass loss driven by stellar feedback may not be 
an effective mechanism to flatten the central cusp. 

Subject headings: galaxies: dwarf — galaxies: halos — galaxies: kinematics and dynamics — galaxies: 
structure — galaxies: evolution 



1. INTRODUCTION 

Cold dark matter (CDM) cosmology is the standard 
paradigm of structure formation in the universe. How- 
ever, it has several serious-unsolved problems. Recent 
observations of nearby dwarf galaxies and low surface 
brightness galaxies (LSBs) have revealed that the den- 
sity profile of the DM halo is constant at the center of 
such galaxies (e.g., Moore 1994; Burkert 1995; de Blok 
et al. 2001; Swaters et al. 2003; Spekkens et al. 2005). 
In contrast, cosmological A^-body simulations based on 
coUisionless CDM have always predicted a steep power- 
law mass-density distribution at the center of CDM ha- 
los (e.g., Navarro et al. 1997; Fukushige & Makino 1997; 
Moore et al. 1998). This discrepancy is the well-known 
core-cusp problem of the CDM scenario. To solve this 
problem, so far, numerical simulations and analytical 
methods used to study the dynamical response of DM ha- 
los to the variance of the gravitational potential. Super- 
nova feedback to the interstellar medium plays a signifi- 
cant role on forming less massive galaxies, because such 
galaxies have a shallower gravitational potential than gi- 
ant galaxies. The effects of stellar feedback differ sub- 
stantially from one galaxy to another depending on the 
gravitational potential (Mori et al. 1999; MacLow & Fer- 
rara 1999). In low-mass dwarf galaxies, supernovae blow 
gas out from the galaxy centers. This mass loss makes the 
gravitational potential around the center of the DM ha- 
los shallower, and the mass-density distribution has a flat 
core (Navarro et al. 1996a; Read & Gilmore 2005). These 
previous works have always assumed that mass loss from 
the center of DM halos occurs instantaneously or within 
a fixed period. However, the timescale of mass loss may 
be also a key factor in determining the dynamical re- 
sponse of DM halos. Furthermore, each dwarf galaxy in 
the local group has a unique star formation history (Ma- 
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teo 1998; Weisz et al. 2011). De Souza et al. (2011) 
analytically demonstrated that the energy gain of par- 
ticles flattens the central cusp; however, they assumed 
that each particle receives an equal amount of energy. 
This assumption must be verified by numerical simula- 
tions with sufficient accuracy. Therefore, in this paper, 
we run coUisionless A^-body simulations that focus on 
the dependence of the dynamical response of DM halos 
on the timescale of gravitational potential variance. In 
other words, we study how the dynamics of DM halos de- 
pend on star formation activities. We demonstrate that 
the timescale of mass loss is one of the important factors 
in determining the dynamical response of DM halos. Fi- 
nally, we discuss the surprising result: the mass loss of 
the baryon is not an effective mechanism for flattening 
the central cusp. 

The structure of this paper is as follows. In §2, we 
describe the numerical simulations. In §3, we show the 
simulation results. Finally, we discuss the results in §4. 

2. NUMERICAL MODELS 

2.1. The DM halo model 

The density distribution of a DM halo obtained from 
A^-body simulations based on CDM cosmology is well 
fitted by the following: 



PDM(f) 



r"(r-f i?DM)3^ 



(1) 



where r is the distance from the center of a DM halo, 
a is a power-low index, and po and i?DM are the scale 
density and scale length of the DM halo, respectively. 
The density distribution changes from p oc r"" in the 
center of a DM halo (r < Rum) to p cx its out- 
skirts (r > Rbm)- Here, a — 1.0 corresponds to 
the Navarro-Frenk- White (NEW) model (Navarro et al. 
1996b; Navarro et al. 1997) and a = 1.5 corresponds to 
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the Fukushige-Makino-Moorc (FMM) model (Fukushige 
& Makino 1997; Moore et al. 1999). To generate equi- 
librium TV-body systems, we use the fitting formulation 
of the distribution function (DF) proposed by Widrow 
(2000). In this case, the DF only depends on energy, and 
the velocity dispersion of the system is isotropic. 

2.2. The baryon model 

Hernquist model (Hernquist 1990) is frequently used 

to describe bulges or elliptical galaxies, because its sur- 
face brightness profile matches dc Vaucouleurs law (de 
Vaucouleurs 1948). Therefore, we adopt the Hernquist 
potential to describe the baryon potential around the 
center of DM halos. The external potential is given by 

where G is the gravitational constant, and Mb and Rh 
are the baryon mass and scale length of the external po- 
tential, respectively. To generate the initial condition, 
we relax the equilibrium A''-body system quoted above 
in the external baryon potential for ~ 30 ta, where td{r) 
is the dynamical time defined by 



TABLE 1 
Summary of simulation run 



37r 



32Gp(< rd) 



(3) 



where p{< r) is the average density of a DM halo within 
a radius r = 0.2 kpc throughout this paper. Then, to 
represent baryon ejection from the center of galaxies, 
we change the baryon mass Mb = Mb,tot(l — V^out), 
where Mb. tot and Tout arc the total baryon mass and the 
timescale of mass loss. For the instantaneous mass-loss 
model, we simply set Mb = 0. 

2.3. Description of runs 

In this paper, we simulate the dynamical response of a 
DM halo with the virial mass Myir = IO^A'/q, the virial 
radius i?vir = 10 kpc, and the scale length Rbm = 2 kpc. 
In this case, the mean dynamical time within 200 pc is 

= 4 Myr for the FMM model and 10 Myr for the NFW 
model. The total baryon mass is Mb, tot = 1-7 x IO^Mq, 
and the scale length of the baryon potential is i?b = 0.04 
kpc. Here, we assume that the total baryon mass in- 
cluded in the dwarf galaxies is initially 17% of the halo 
mass. This fraction matches the cosmic baryon frac- 
tion /b obtained from Wilkinson Microwave Anisotropy 
Probe (WMAP) observations (Spergel et al. 2007; Ko- 
matsu et al. 2009). The adopted parameters for each of 
our simulations are given in Table 1. Throughout this 
paper, the opening parameter of the Barnes-Hut tree al- 
gorithm (Barnes & Hut 1986) is9 = 0.8 and the softening 
parameter is e = 0.008 kpc. 

3. SIMULATION RESULTS 

In this section, we present the simulation results show- 
ing the evolutionary processes in the instantaneous and 
quasi-adiabatic mass-loss models. 

3.1. Dynamical evolution of DM halos 

We show how the density profile of a DM halo evolves 
after mass loss for the instantaneous mass-loss case with 
the FMM initial condition (MLl) at the left panels in 
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UPl 


FMM 


1,048,576 


495 td 






UP2 


FMM 


16,384 


7.74 td 






UP3 


NFW 


1,048,576 


148 td 






UP4 


NFW 


16,384 


2.32 id 






MLl 


FMM 


1,048,576 


495 id 


instantaneous 


0.85 


ML2 


FMM 


1,048,576 


495 td 


1 td 


1.20 


ML3 


FMM 


1,048,576 


495 id 


10 id 


1.44 


ML4 


FMM 


1,048,576 


495 td 


50 td 


1.46 


ML5 


NFW 


1,048,576 


148 id 


instantaneous 


0.42 


ML6 


NFW 


1,048,576 


148 td 


50 td 


0.89 


Note. 


— ML indicates mass- 


loss runs 


, and UP indicates unper- 



turbed runs to examine the stability of the A''-body system and the 
effects of two-body relaxation. T^ei is a two-body relaxation time 
within 0.2 kpc. 



Fig. 1. During the first ~ 30 td (see Fig. lb), the sys- 
tem expands and the central cusp is flattened to the core. 

In this phase, the bulk of particles in the center of DM 
halo moves outward, and then the self-gravity of the DM 
halo slows the particles. After 30 t^, the particles lose 
their radial velocity and fall back into the center. In Fig. 
Ic-ld, these particles assemble at the center of the halo, 
and the cusp is regenerated. This cusp regenerating pro- 
cess occurs from the center to the outward, and finally, 
the system reaches the quasi-cquilibrium state (see Fig 
le) . We emphasize that the central cusp has been regen- 
erated, but the power-law index of density at the center 
a = 0.85 is smaller than the initial condition a = 1.5. A 
similar phenomenon occurs for the instantaneous mass- 
loss case with the NFW initial condition (MLS): the den- 
sity at the center is a = 0.42 for the initial condition 
a = 1.0. Therefore, instantanciOTis mass- loss flattens the 
central cusp in the density distribution to some degree. 

In the right panels of Fig. 1, we show how the density 
profile of a DM halo evolves during and after mass loss for 
the slow (quasi-adiabatic) mass-loss case with the FMM 
initial condition (ML4). Before the mass- loss operation, 
the DM halo has a high central density because of the 
baryon potential; this density decreases with decreasing 
baryon mass. Slow mass-loss does not significantly flat- 
ten the central cusp in the density distribution. More- 
over, the system always keeps the virial state for the slow 
mass- loss process, and the final density profile of the DM 
halo recovers the FMM initial condition. A similar phe- 
nomenon occurs for the adiabatic mass-loss case with the 
NFW initial condition (ML6). 

3.2. Final states of DM halos 

Figure 2 shows the results of simulations for the FMM 

model (left panels) and the NFW model (right panels). 
The upper left (right) panel shows the density profile of 
a DM halo for the unperturbed model with the FMM 
(NFW) initial condition, respectively. The thick-solid 
line corresponds to the high-resolution run (UPl and 
UP3) and the thin-solid line corresponds to the low- 
resolution run (UP2 and UP4). In unperturbed runs with 
only a small number of particles, even though no exter- 
nal potential exists, a cusp-to-core transition caused by 
two-body relaxation occurs. Previous studies could not 
rule out this effect, because they did not use a sufficient 
number of particles in the central region. However, two of 
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Fig. 1. — Evolution of the density profiles of a DM halo for the 
instantaneous mass-loss model (MLl; left panels) and the adiabatic 
mass-loss model (ML4; right panels). The top panel (a) shows 
the density profile of the quasi-cquilibrium state after adding the 
external potential on the initial FMM model. The other panels 
show the density profiles of a DM halo at 15 (b and f), 30 <d (c 
and g), 50 (d and h) and 110 (e and i), respectively. Dashed- 
line represents the FMM initial condition. 

unperturbed runs clearly show that we have a sufficient 
number of particles to suppress the effect of two-body 
relaxation and the system remains stable at least several 
100 dynamical times. 

The middle panels in Fig. 2 show that the density pro- 
file of a DM halo has reached the quasi-cquilibrium state 
after mass loss in ML runs. The thick-solid lines repre- 
sent the instantaneous mass-loss model (MLl and MLS) 
and the thin-solid line corresponds Tout = 50td (AfL4 
and ML6). Wc fit the quasi-equilibrium density profiles 
for all of ML runs using the least-squares method, and 
wc find the following values for the parameter a. For the 
FMM model, a is 0.85 for the instantaneous mass-loss 
(MLl), 1.20 for Tout = id (ML2). 1.44 for Tout = 10 ta 
(ML3), and 1.46 for Tout = 50 t,i (ML4). In Table 1, we 
summarize the resultant a including for the NFW model. 
It is clear that mass loss occurs in a short timescale; in 
other words, intense stellar activities cause the density 
profiles of DM halos to become shallower. In the slow 
mass-loss cases (MLS, ML4 and ML6; Tout ^ 10 td), a is 
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Fig. 2. — The upper left jjancl shows the density profile of a 
DM halo for the unperturbed (UP) model (solid lines) with the 
FMM initial condition (long-dashed line). The thick-solid line rep- 
resents A*" = 1,048,576 (UPl) and the thin-solid line represents 
N = 16, 384 (UP2). The middle left (right) panel shows the quasi- 
equilibrium density profile of a DM halo for the FMM (NFW) 
model with the mass-loss. The thick-solid line represents the in- 
stantaneous mass-loss model (MLl and ML5), and the thin-solid 
line corresponds Tout = SOtj (ML4 and ML6). The lower left 
(right) panel shows the variance of the total energy of the system 
for the FMM model (NFW model). Thick-solid line corresponds 
the instantaneous mass-loss model, and thin-solid line corresponds 
Tout = 50 td ) where Eq and E are the total energies of DM halos 
at the initial and given times, respectively. 

approximately conserved. Therefore, to flatten the cen- 
tral cusp needs intense stellar activity; however, even for 
cases MLl and ML5 that have such intense stellar activ- 
ity, a is still larger than the typical observational values 
of a ~ 0.2 - 0.3 (Spekkens et al. 2005). 

The lower panels of Fig. 2 show the amount of en- 
ergy gained from the external baryon potential. While 
the baryon potential is added into DM halos, particles in 
the system acquire kinetic energy to balance the potc!n- 
tial energy; conversely, during mass loss, particles lose 
kinetic energy. After mass-loss operation has finished, 
the total energy conserves. The lower panels of Fig. 1 
demonstrate that after mass loss in the instantaneous 
mass-loss cases (MLl and ML5; thick solid lines), DM 
halos gain more energy than that gained in the slower 
mass-loss cases (ML2 and ML6: thin solid lines). There- 
fore, the shorter the timescale of mass loss, the farther 
the DM halos expand and the flatter they become. 

4. DISCUSSION 

Previous numerical attempts to solve the core-cusp 
problem of the CDM scenario (Navarro et al. 1996a; 
Read & Gilmore 2005) have been limited by two factors: 

(1) they did not consider the timescale of mass loss and 

(2) they could not remove the effects of the two-body re- 
laxation, because they did not use the sufficient number 
of particles in the center of the DM halos. Our current 
numerical study overcomes these limitations by consid- 
ering the timescale of mass loss and including a larger 
number of particles in the center of DM halos. We flnd 
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Fig. 3. — The distribution function of DM hales on the nor- 
mahzed particle energy. This energy is normalized by multiplying 
by the normalization factor A. Shown are the initial conditions 
for the FMM model (thick dashed hne; A = -0.70), the quasi- 
equilibrium state before the mass-loss operation (thin dashed line; 
A = —0.095), the quasi-equilibrium state for the instantaneous 
mass-loss model MLl (thick solid line; A = —2.0) and the quasi- 
adiabatic mass-loss model ML4 (thin solid line; A = —1.05). 

that the timescale of mass loss is an important factor 
affecting DM halo dynamics. 

Our study determines the dynamical response of DM 
halos to changes in the gravitational potential induced by 
stellar activities that remove the baryons. We especially 
focus on the timescale of gravitational variance. We find 
that the central cusp of the DM halo is flatter when mass 
loss occurs over a short time than when it occurs over a 
long time; this result is consistent with the analytical 
insight of Hills (1980). However, the power-law index of 
the central cusp a is still larger than typical observational 
values; in other words, the central cusp remains even 
for instantaneous mass-loss models. Therefore, mass loss 
may not be an effective mechanism to flatten the central 
cusp, at least in spherical systems. This suggests that 
DM halos may be stable to mass loss. 

Figure 3 shows the DF for the FMM initial condition, 
the quasi-equilibrium state with the external baryon po- 
tential, MLl and ML4. The shape of the DFs after mass 
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loss is similar to the initial condition, suggesting that DM 
halos return to the initial condition not only in spatial- 
space but also in phase-space; therefore, DM halos are 
stable to mass loss. 

In this study, we find that a temporal cusp-to-core 
transition occurs for instantaneous mass- loss. After a 
few dozen dynamical times, the cusp regenerates from the 
center. However, this cusp has a smaller a than the ini- 
tial conditions, suggesting that the cusp has been some- 
what flattened by mass loss. In addition, we find that 
not only the density profile but also the energy distribu- 
tion almost returns to the initial condition for adiabatic 
mass-loss. The physical mechanism of this cusp regen- 
eration has not been revealed, but it may aid in under- 
standing the formation of the universal density profiles of 
CDM halos. Only in the linear perturbation regime, the 
Doremus-Feix-Baumann theorem (Doremus et al. 1971) 
states that coUisionless TV-body systems with DF f{E), 
df{E)/dE < are stable to radial perturbation. How- 
ever, it is open question in the non-linear regime. In 
the series of subsequent studies, we will elucidate the 
detained mechanism of the recovering cusp, and exam- 
ine more realistic model including non-spherical analysis 
and gasdynamics. 
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